ggvote_plot_mfx2 <- function(ames_data,
                             desc_data,
                             model_ids,
                             add_descriptives = TRUE,
                             d = "cmr_arm_cwu",
                             x1  = "Homeownership",
                             x2,
                             x2_continuous = TRUE,
                             subset_x1 = NULL,
                             ylim = rep(NA_real_, 2L),
                             ytitle = "Marginal Effect",
                             xtitle = "",
                             font_size = 15L) {
  ## Data
  if (x2_continuous) {
    data <- ames_data %>%
      dplyr::filter(model_id == model_ids[1]) %>%
      dplyr::bind_rows(ames_data %>%
                         dplyr::filter(model_id == model_ids[2])) %>%
      dplyr::filter(factor == d) %>%
      dplyr::mutate(
        Homeownership = ifelse(model_id == model_ids[1],
                               "Renters",
                               "Owners") %>%
          factor(.,
                 levels = c("Renters", "Owners"))
      )
    # x2_breaks_1 <-
    #   sort(unique(data[data$model_id == model_ids[1], x2]))
    # x2_breaks_2 <-
    #   sort(unique(data[data$model_id == model_ids[2], x2]))
    # x2_breaks_1_mod <- c(x2_breaks_1[1] - abs(diff(x2_breaks_1[1:2])),
    #                      x2_breaks_1,
    #                      x2_breaks_1[length(x2_breaks_1)] +
    #                        abs(diff(x2_breaks_1[(length(x2_breaks_1) - 1):length(x2_breaks_1)])))
    # x2_breaks_2_mod <- c(x2_breaks_2[1] - abs(diff(x2_breaks_2[1:2])),
    #                      x2_breaks_2,
    #                      x2_breaks_2[length(x2_breaks_2)] +
    #                        abs(diff(x2_breaks_2[(length(x2_breaks_2) - 1):length(x2_breaks_2)])))
    desc <- desc_data %>%
      dplyr::filter(model_id %in% model_ids) %>%
      dplyr::filter(variable == x2) %>%
      # dplyr::mutate(!!as.name(x2) := ifelse(
      #   model_id == model_ids[1],
      #   cut(breaks, breaks = c(-Inf, x2_breaks_1, Inf)),
      #   cut(breaks, breaks = c(-Inf, x2_breaks_2, Inf))
      # )) %>%
      group_by(model_id) %>%
      dplyr::mutate(
        x2_lower = ifelse(row_number() < n(),
                          breaks,
                          NA_real_),
        x2_upper = ifelse(row_number() < n(),
                          dplyr::lead(breaks),
                          NA_real_)
      ) %>%
      dplyr::filter(not(is.na(x2_upper))) %>%
      dplyr::mutate(
        Homeownership = ifelse(model_id == model_ids[1],
                               "Renters",
                               "Owners") %>%
          factor(.,
                 levels = c("Renters", "Owners"))
      )
  } else {
    join_vec <- c("model_id", "levels")
    names(join_vec) <- c("model_id", x2)
    data <- ames_data %>%
      dplyr::filter(model_id == model_ids[1]) %>%
      dplyr::bind_rows(ames_data %>%
                         dplyr::filter(model_id == model_ids[2])) %>%
      dplyr::filter(factor == d) %>%
      dplyr::left_join(
        desc_data %>%
          dplyr::filter(model_id %in% model_ids) %>%
          dplyr::filter(variable %in% x2) %>%
          dplyr::group_by(model_id, levels) %>%
          dplyr::summarize(counts = sum(counts, na.rm = TRUE)),
        by = join_vec
      ) %>%
      dplyr::mutate(
        Homeownership = ifelse(model_id == model_ids[1],
                               "Renters",
                               "Owners") %>%
          factor(.,
                 levels = c("Renters", "Owners"))
      ) %>%
      dplyr::ungroup()
  }
  
  ## Subset by x1
  if (not(is.null(subset_x1))) {
    if (all(is.na(ylim))) {
      ylim <- c(min(data$lower, na.rm = TRUE),
                max(data$upper, na.rm = TRUE))
    }
    
    data <- data %>%
      filter(!!as.name(x1) == subset_x1)
    
    if (x2_continuous) {
      desc <- desc %>%
        filter(!!as.name(x1) == subset_x1)
    }
  }
  
  ## Marginal effect
  vote_plot <- data %>%
    ggplot(aes(x = !!as.name(x2),
               y = AME)) +
    geom_point(position = position_dodge(width = 1)) +
    geom_errorbar(aes(ymin = lower,
                      ymax = upper),
                  position = position_dodge(width = 1),
                  width = 0) +
    geom_hline(yintercept = 0, col = "red") +
    facet_grid(. ~ get(eval(x1))) +
    theme(panel.spacing = unit(0.75, "lines"))  +
    theme(text = element_text(size = font_size)) +
    ylim(ylim) +
    ylab(ytitle) +
    xlab(xtitle)
  
  ## Descriptives
  if (add_descriptives) {
    if (is.null(subset_x1)) {
      effective_ymin <-
        ggplot_build(vote_plot)$layout$panel_scales_y[[1]]$range$range[1]
      effective_ymax <-
        ggplot_build(vote_plot)$layout$panel_scales_y[[1]]$range$range[2]
    } else {
      effective_ymin <- ylim[1]
      effective_ymax <- ylim[2]
    }
    
    if (x2_continuous) {
      vote_plot <- vote_plot +
        geom_rect(
          data = desc %>%
            dplyr::rowwise() %>%
            dplyr::mutate(x2_mid = mean(c(
              x2_lower, x2_upper
            ))) %>%
            dplyr::ungroup() %>%
            dplyr::mutate(
              counts = dplyr::coalesce(counts, 0L),!!as.name(x2) := NA_real_,
              AME = NA_real_
            ),
          aes(
            xmin = x2_lower,
            xmax = x2_upper,
            ymin = effective_ymin,
            ymax = ((counts / max(counts)) *
                      abs(effective_ymax - effective_ymin)) +
              effective_ymin
          ),
          alpha = .15,
          color = adjustcolor("gray20", alpha.f = .25)
        ) +
        xlim(range(data[[x2]]))
    } else {
      vote_plot <- vote_plot +
        geom_segment(
          aes(
            x = !!as.name(x2),
            xend = !!as.name(x2),
            y = effective_ymin,
            yend = ((counts / max(counts)) *
                      abs(effective_ymax - effective_ymin)) +
              effective_ymin
          ),
          size = 25,
          alpha = .15
        )
    }
  }
  
  ## Return
  return(vote_plot)
}
